!>\file sgscloud_radpre.F90
!!  Contains the preliminary (interstitial) work to the call to the radiation schemes:
!!    1) Backs up the original qc & qi
!!    2) Adds the partioning of convective condensate into liqice/ice for effective radii
!!    3) Adds the subgrid clouds mixing ratio and cloud fraction to the original (resolved-
!!       scale) qc, qi and cloud fraction coming from the microphysics scheme.
!!    4) Recompute the diagnostic high, mid, low, total and bl clouds to be consistent with radiation

      module sgscloud_radpre

      contains

!> \defgroup sgsradpre_group sgscloud_radpre_run Module
!! This interstitial code adds the subgrid clouds to the resolved-scale clouds
!! if there is no resolved-scale clouds in that particular grid box. It can also
!! specify a cloud fraction for resolved-scale clouds as is done currently when
!! using MYNN-EDMF. For clouds coming from the convection schemes (in this case
!! only used by GF scheme), two cloud fraction options are available:
!! Xu-Randall (XR1996) or Chaboureau and Bechtold (CB2005), chosen by the
!! switch "conv_cf_opt" = 0: CB2005, 1: XR1996.
!!
!> \section arg_table_sgscloud_radpre_run Argument Table
!! \htmlinclude sgscloud_radpre_run.html
!!
!!    cloud array description:                                          !
!!          clouds(:,:,1)  -  layer total cloud fraction                !
!!          clouds(:,:,2)  -  layer cloud liq water path                !
!!          clouds(:,:,3)  -  mean effective radius for liquid cloud    !
!!          clouds(:,:,4)  -  layer cloud ice water path                !
!!          clouds(:,:,5)  -  mean effective radius for ice cloud       !
!!          clouds(:,:,6)  -  layer rain drop water path                !
!!          clouds(:,:,7)  -  mean effective radius for rain drop       !
!!          clouds(:,:,8)  -  layer snow flake water path               !
!!          clouds(:,:,9)  -  mean effective radius for snow flake
!!
!>\section sgscloud_radpre_mod  SGS Cloud Scheme Pre General Algorithm
      subroutine sgscloud_radpre_run(    &
           im,dt,levs,                   &
           flag_init,flag_restart,       &
           con_g, con_pi, eps, epsm1,    &
           r_v, cpv, rcp,                &
           xlv, xlf, cp,                 &
           do_mynnedmf,                  &
           qc, qi, qv, T3D, P3D, exner,  &
           qr, qs, qg,                   &
           qci_conv,ud_mf,               &
           imfdeepcnv, imfdeepcnv_gf,    &
           qc_save, qi_save, qs_save,    &
           qc_bl,qi_bl,cldfra_bl,        &
           delp,clouds1,clouds2,clouds3, &
           clouds4,clouds5,              &
           clouds8,clouds9,slmsk,        &
           nlay, plyr, xlat, dz,de_lgth, &
           cldsa,mtopa,mbota,            &
           imp_physics, imp_physics_gfdl,&
           iovr,                         &
           errmsg, errflg                )

      use machine , only : kind_phys
      use module_radiation_clouds, only : gethml
      use radcons,  only: qmin          ! Minimum values for various calculations
      use funcphys, only: fpvs          ! Function to compute sat. vapor pressure over liq.
!------------------------------------------------------------------- 
      implicit none
!------------------------------------------------------------------- 
      ! Interface variables
      real(kind=kind_phys), intent(in) :: con_g, con_pi, eps, epsm1
      real(kind=kind_phys), intent(in) :: r_v, cpv, rcp
      real(kind=kind_phys), intent(in) :: xlv, xlf, cp
      real(kind=kind_phys), intent(in) :: dt
      real :: xls, xlvcp, xlscp !derived below
      real(kind=kind_phys)             :: gfac
      integer,             intent(in)  :: im, levs, imfdeepcnv, imfdeepcnv_gf, &
           &               nlay, imp_physics, imp_physics_gfdl
      logical,             intent(in)  :: flag_init, flag_restart, do_mynnedmf

      real(kind=kind_phys), dimension(:,:), intent(inout) :: qc, qi
      real(kind=kind_phys), dimension(:,:), intent(inout) :: qr, qs, qg
      ! qci_conv only allocated if GF is used
      real(kind=kind_phys), dimension(:,:), intent(inout) :: qci_conv
      real(kind=kind_phys), dimension(:,:), intent(in)    :: ud_mf
      real(kind=kind_phys), dimension(:,:), intent(in)    :: T3D,delp
      real(kind=kind_phys), dimension(:,:), intent(in)    :: qv,P3D,exner
      real(kind=kind_phys), dimension(:,:), intent(inout) ::  &
           &         clouds1,clouds2,clouds3,clouds4,clouds5, &
           &         clouds8,clouds9
      real(kind=kind_phys), dimension(:,:), intent(inout) :: qc_save, qi_save, qs_save
      real(kind=kind_phys), dimension(:,:), intent(in)    :: qc_bl, qi_bl, cldfra_bl
      real(kind=kind_phys), dimension(:),   intent(in)    :: slmsk, xlat, de_lgth
      real(kind=kind_phys), dimension(:,:), intent(in)    :: plyr, dz      
      real(kind=kind_phys), dimension(:,:), intent(inout) :: cldsa
      integer,              dimension(:,:), intent(inout) :: mbota, mtopa
      integer,                              intent(in)    :: iovr
      character(len=*), intent(out) :: errmsg
      integer,          intent(out) :: errflg
      ! Local variables
      ! pressure limits of cloud domain interfaces (low,mid,high) in mb (0.1kPa)
      real(kind=kind_phys) :: ptop1(im,3+1)  !< pressure limits of cloud domain interfaces
      real(kind=kind_phys) :: ptopc(3+1,2 )  !< pressure limits of cloud domain interfaces
                                             !! (low, mid, high) in mb (0.1kPa)
      data ptopc / 1050., 650., 400., 0.0,  1050., 750., 500., 0.0 /
      real(kind=kind_phys), dimension(im,nlay) :: cldcnv
      real(kind=kind_phys), dimension(im)      :: rxlat
      real(kind=kind_phys) :: Tc, Tk, liqfrac, iwc, ice_frac, snow_frac
      integer              :: i, k, id
      ! DH* 20200723 - see comment at the end of this routine around 'gethml'
      real(kind=kind_phys), dimension(im,nlay) :: alpha_dummy
      ! *DH

      ! PARAMETERS FOR RANDALL AND XU (1996) CLOUD FRACTION
      real, parameter  ::  coef_p = 0.25, coef_gamm = 0.49, coef_alph = 100.
      real :: rhgrid,h2oliq,qsat,tem1,tem2,clwt,es,onemrh,value

      !Chaboureau and Bechtold (2002 and 2005)
      real :: a, f, sigq, qmq, qt, xl, tlk, th, thl, rsl, cpm, cb_cf

      !Option to convective cloud fraction
      integer, parameter :: conv_cf_opt = 0  !0: C-B, 1: X-R

      ! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

      ! some derived variables from incoming constants:
      xls=xlv+xlf
      xlvcp=xlv/cp
      xlscp=(xlv+xlf)/cp

      !write(0,*)"=============================================="
      !write(0,*)"in SGSCLoud_RadPre"
      gfac=1.0e5/con_g
      if (flag_init .and. (.not. flag_restart)) then
        !write (0,*) 'Skip this flag_init = ', flag_init
        ! return
        ! Need default cloud fraction when MYNN is not used: Resort to
        ! Xu-Randall (1996).
        ! cloud fraction =
        ! {1-exp[-100.0*qc/((1-RH)*qsat)**0.49]}*RH**0.25
        do k = 1, levs
          do i = 1, im
            if ( qi(i,k) > 1E-7 .OR. qc(i,k) > 1E-7 ) then
              es     = min( p3d(i,k),  fpvs( t3d(i,k) ) )  ! fpvs and prsl in pa
              qsat   = max( QMIN, eps * es / (p3d(i,k) + epsm1*es) )
              rhgrid = max( 0., min( 1., qv(i,k)/qsat ) )
              h2oliq = qc(i,k) + qi(i,k) + qr(i,k) + qs(i,k) + qg(i,k)   ! g/kg
              clwt   = 1.0e-6 * (p3d(i,k)*0.00001)

              if (h2oliq > clwt) then
                onemrh= max( 1.e-10, 1.0-rhgrid )
                tem1  = min(max((onemrh*qsat)**0.49,0.0001),1.0)  !jhan
                tem1  = 100.0 / tem1
                value = max( min( tem1*(h2oliq-clwt), 50.0 ), 0.0 )
                tem2  = sqrt( sqrt(rhgrid) )

                clouds1(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
              endif

            endif
          enddo
        enddo

      else ! timestep > 1 or restart

        ! Back-up microphysics cloud information:
        do k = 1, levs
          do i = 1, im
            qc_save(i,k) = qc(i,k)
            qi_save(i,k) = qi(i,k)
            qs_save(i,k) = qs(i,k)
          end do
        end do

        if ( do_mynnedmf ) then

          ! add boundary layer clouds - Note: now the temperature-dependent sorting of
          ! ice and water subgrid-scale clouds is done inside the MYNN-EDMF

          do k = 1, levs
            do i = 1, im

              !if (imp_physics == imp_physics_gfdl) then
              !  ! only complement the GFDL cloud fractions
              !  if (clouds1(i,k) < 0.01 .and. cldfra_bl(i,k) > 0.01) then
              !    clouds1(i,k) = cldfra_bl(i,k)
              !  endif
              !else
                 clouds1(i,k) = cldfra_bl(i,k)
              !endif

              if (qc(i,k) < 1.e-6 .and. cldfra_bl(i,k)>0.001) then
                qc(i,k) = qc_bl(i,k)*cldfra_bl(i,k)

                !eff radius cloud water (microns) from Miles et al. (2007)
                if (nint(slmsk(i)) == 1) then !land
                  if(qc(i,k)>1.E-8)clouds3(i,k)=5.4
                else
                  if(qc(i,k)>1.E-8)clouds3(i,k)=9.6
                endif

                !calculate the liquid water path using additional BL clouds
                clouds2(i,k) = max(0.0, qc(i,k) * gfac * delp(i,k))
              endif

              Tc = T3D(i,k) - 273.15
              !crudely split frozen species into 50% ice and 50% snow below
              !~700 mb and decrease snow to zero by ~300 mb 
              snow_frac = min(0.5, max((p3d(i,k)-30000.0),0.0)/140000.0)
              ice_frac  = 1.0 - snow_frac
              if (qi(i,k) < 1.e-8 .and. cldfra_bl(i,k)>0.001) then
                qi(i,k) = ice_frac*qi_bl(i,k)*cldfra_bl(i,k)

                !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 6b)
                if(qi(i,k)>1.E-8)clouds5(i,k)=max(173.45 + 2.14*Tc, 20.)
                !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 8b)
                !iwc = qi(i,k)*1.0e6*rho(i,k)
                !IF(qi(i,k)>1.E-8)clouds5(i,k)=MAX(139.7 + 1.76*Tc + 13.49*LOG(iwc), 20.)

                !calculate the ice water path using additional BL clouds
                clouds4(i,k) = max(0.0, qi(i,k) * gfac * delp(i,k))
              endif

              if (qs(i,k) < 1.e-8 .and. cldfra_bl(i,k)>0.001) then
                qs(i,k) = snow_frac*qi_bl(i,k)*cldfra_bl(i,k)

                !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 6b)
                if(qs(i,k)>1.E-8)clouds9(i,k)=max(2.*(173.45 + 2.14*Tc), 50.)

                !calculate the snow water path using additional BL clouds
                clouds8(i,k) = max(0.0, qs(i,k) * gfac * delp(i,k))
              endif
            enddo
          enddo

        elseif (imp_physics /= imp_physics_gfdl) then 

          ! Non-MYNN cloud fraction AND non-GFDL microphysics, since both
          ! have their own cloud fractions. In this case, we resort to
          ! Xu-Randall (1996).
          ! cloud fraction =
          ! {1-exp[-100.0*qc/((1-RH)*qsat)**0.49]}*RH**0.25
          do k = 1, levs
            do i = 1, im
              if ( qi(i,k) > 1E-7 .OR. qc(i,k) > 1E-7 ) then

                es     = min( p3d(i,k),  fpvs( t3d(i,k) ) )  ! fpvs and prsl in pa
                qsat   = max( QMIN, eps * es / (p3d(i,k) + epsm1*es) )
                rhgrid = max( 0., min( 1., qv(i,k)/qsat ) )
                h2oliq = qc(i,k) + qi(i,k) + qr(i,k) + qs(i,k) + qg(i,k) ! g/kg
                clwt   = 1.0e-6 * (p3d(i,k)*0.00001)

                if (h2oliq > clwt) then
                  onemrh= max( 1.e-10, 1.0-rhgrid )
                  tem1  = min(max((onemrh*qsat)**0.49,0.0001),1.0)  !jhan
                  tem1  = 100.0 / tem1
                  value = max( min( tem1*(h2oliq-clwt), 50.0 ), 0.0 )
                  tem2  = sqrt( sqrt(rhgrid) )

                  clouds1(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
                endif

              endif
            enddo
          enddo

        endif ! end MYNN or OTHER choice for background clouds fractions

        ! At this point, we have cloud properties for all non-deep convective clouds.
        ! So now we add the convective clouds:

        if (imfdeepcnv == imfdeepcnv_gf) then
          do k = 1, levs
            do i = 1, im
              !if ( qci_conv(i,k) > 0. .AND. (qi(i,k) < 1E-7 .AND. qc(i,k) < 1E-7 ) ) then
              if ( qci_conv(i,k) > 0. ) then
                Tk = T3D(i,k)
                Tc = Tk - 273.15

                !Partition the convective clouds into water & frozen species
                liqfrac = min(1., max(0., (Tk-244.)/29.))
                qc(i,k) = qc(i,k)+qci_conv(i,k)*liqfrac
                !split ice & snow 50-50%
                qi(i,k) = qi(i,k)+0.5*qci_conv(i,k)*(1. - liqfrac)
                qs(i,k) = qs(i,k)+0.5*qci_conv(i,k)*(1. - liqfrac)

                !eff radius cloud water (microns)
                if (nint(slmsk(i)) == 1) then !land
                  if(qc(i,k)>1.E-8)clouds3(i,k)=5.4
                else
                  !from Miles et al.
                  if(qc(i,k)>1.E-8)clouds3(i,k)=9.6
                endif
                !from Mishra et al. (2014, JGR Atmos), assume R_sno = 2*R_ice
                if(qi(i,k)>1.e-8)clouds5(i,k)=max(     173.45 + 2.14*Tc , 20.)
                if(qs(i,k)>1.e-8)clouds9(i,k)=max(2.0*(173.45 + 2.14*Tc), 50.)

                if ( conv_cf_opt .eq. 0 ) then
                   !print *,'Chab-Bechtold cloud fraction used'
                   !  clouds1(i,k) = cldfra_bl(i,k)

                   !Alternatively, use Chaboureau-Bechtold (CB) convective component
                   !Based on both CB2002 and CB2005.
                   xl  = xlv*liqfrac + xls*(1.-liqfrac)  ! blended heat capacity
                   tlk = t3d(i,k) - xlvcp/exner(i,k)*qc(i,k) &
                       &          - xlscp/exner(i,k)*qi(i,k)! liquid temp
                   ! get saturation water vapor mixing ratio at tl and p
                   es  = min( p3d(i,k), fpvs( tlk ) )   ! fpvs and prsl in pa
                   qsat= max( QMIN, eps*es / (p3d(i,k) + epsm1*es) )
                   rsl = xl*qsat / (r_v*tlk**2)   ! slope of C-C curve at t = tl 
                                                  ! CB02, Eqn. 4
                   qt  = qc(i,k) + qi(i,k) + qv(i,k) !total water
                   cpm = cp + qt*cpv              ! CB02, sec. 2, para. 1
                   a   = 1./(1. + xl*rsl/cpm)     ! CB02 variable "a"
                   !Now calculate convective component of the cloud fraction:
                   if (a > 0.0) then
                      f = min(1.0/a, 4.0)         ! f is the vertical profile
                   else                           ! scaling function (CB2005)
                      f = 1.0
                   endif
                   sigq = 1.5E-3 * ud_mf(i,k)/dt * f
                   !sigq = 3.E-3 * ud_mf(i,k)/dt * f
                   sigq = SQRT(sigq**2 + 1e-10)   ! combined conv + background components
                   qmq  = a * (qt - qsat)         ! saturation deficit/excess;
                                                  !   the numerator of Q1
                   cb_cf= min(max(0.5 + 0.36 * atan(1.55*(qmq/sigq)),0.01),0.99)
                   if (do_mynnedmf .and. qmq .ge. 0.0) then
                      ! leverage C-B stratus clouds from MYNN in saturated conditions
                      clouds1(i,k) = 0.5*(clouds1(i,k) + cb_cf)
                   else                           ! unsaturated
                      clouds1(i,k) = cb_cf
                   endif
                else
                   !print *,'GF with Xu-Randall cloud fraction'
                   ! Xu-Randall (1996) cloud fraction
                   es     = min( p3d(i,k),  fpvs( t3d(i,k) ) )  ! fpvs and prsl in pa
                   qsat   = max( QMIN, eps*es / (p3d(i,k) + epsm1*es) )
                   rhgrid = max( 0., min( 1.00, qv(i,k)/qsat ) )
                   h2oliq = qc(i,k) + qi(i,k) + qr(i,k) + qs(i,k) + qg(i,k) ! g/kg
                   clwt   = 1.0e-6 * (p3d(i,k)*0.00001)

                   if (h2oliq > clwt) then
                      onemrh= max( 1.e-10, 1.0-rhgrid )
                      tem1  = min(max((onemrh*qsat)**0.49,0.0001),1.0)  !jhan
                      tem1  = 100.0 / tem1
                      value = max( min( tem1*(h2oliq-clwt), 50.0 ), 0.0 )
                      tem2  = sqrt( sqrt(rhgrid) )

                      clouds1(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
                   else
                      clouds1(i,k) = 0.0
                   endif
                   !print*,"XuRandla- cf:",clouds1(i,k)," rh:",rhgrid," qt:",h2oliq
                   !print*,"XuRandlb- clwt:",clwt," qsat:",qsat," p:",p3d(i,k)
                endif ! end convective cf choice
              endif ! qci_conv
            enddo
          enddo
        endif ! imfdeepcnv_gf

      endif ! timestep > 1

      end subroutine sgscloud_radpre_run
      end module sgscloud_radpre
